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We address the phenomenon of drag reduction by dilute polymeric additive to turbulent flows, 
using Direct Numerical Simulations (DNS) of the FENE-P model of viscoelastic flows. It had been 
amply demonstrated that these model equations reproduce the phenomenon, but the results of DNS 
were not analyzed so far with the goal of interpreting the phenomenon. In order to construct a 
useful framework for the understanding of drag reduction we initiate in this paper an investigation 
of the most important modes that are sustained in the viscoelastic and Newtonian turbulent flows 
respectively. The modes are obtained empirically using the Karhunen-Loeve decomposition, allowing 
us to compare the most energetic modes in the viscoelastic and Newtonian flows. The main finding 
of the present study is that the spatial profile of the most energetic modes is hardly changed between 
the two flows. What changes is the energy associated with these modes, and their relative ordering 
in the decreasing order from the most energetic to the least. Modes that are highly excited in one 
flow can be strongly suppressed in the other, and vice versa. This dramatic energy redistribution 
is an important clue to the mechanism of drag reduction as is proposed in this paper. In particular 
there is an enhancement of the energy containing modes in the viscoelastic flow compared to the 
Newtonian one; drag reduction is seen in the energy containing modes rather than the dissipative 
modes as proposed in some previous theories. 



I. INTRODUCTION 

"Drag reduction" refers to the interesting observation 
that the addition of a few tens of parts per million (by 
weight) of long-chain polymers to turbulent fluids can 
bring about a reduction of the friction drag by up to 
80% 0. Obviously, the phenomenon has far reaching 
practical implications besides being challenging from the 
fundamental point of view. In spite of intense interest 
for an extended period of time |2|, ||, ^| , Sreenivasan and 
White recently concluded that "it is fair to say that 
the extensive - and continuing - activity has not produced 
a firm grasp of the mechanisms of drag reduction" . In 
this paper we want to advance on the basis of recent 
Direct Numerical Simulation of model viscoelastic hy- 
drodynamic equations (|, ||, @, §|. Such DNS show un- 
equivocally that drag reduction is reproduced by model 
equations like the FENE-P model. From the theoreti- 
cal point of view this is significant, since it indicates that 
the phenomenon is included in the solutions of the model 
equations. Understanding drag reduction then becomes 
a usual challenge of theoretical physics. 

Our thinking here is motivated in part by a recent anal- 
ysis of the stability of laminar channel flows subject to 
space dependent effective viscosity (9[ [lOj. It turned out 
the even small viscosity gradients can lead to a giant sta- 
bilization of the most unstable modes, both for primary 
and secondary instabilities. In these cases one can un- 
derstand the phenomenon completely by examining the 
energy budget of the putative unstable modes and their 
interaction with the mean flow; the most important ob- 
servation had been that it is the the existence of viscosity 
gradients positioned at a strategic distance from the wall 
which is crucial for the existence of a large effect. It seems 



desirable to do something similar for the viscoelastic tur- 
bulent flows as well (in which the space dependent effec- 
tive viscosity arises due to differential stretching of the 
polymers). But alas, in distinction from primary and sec- 
ondary instabilities, where it is obvious which are the rel- 
evant modes, for the turbulent flow these are not known 
apriori. We therefore decided to first initiate a system- 
atic study of the empirical modes that are sustained in 
the turbulent flow, and then to discuss their interaction 
with the mean flow and with the polymeric additive, their 
stabilization or destabilization when we compare the vis- 
coelastic to the Newtonian flow, and their energy budget. 
In this paper we present the first results of this study. 

We will demonstrate that we can determine with rea- 
sonable accuracy at least the first thirty most energetic 
modes that are sustained in the turbulent flow, for both 
the FENE-P and the Navier-Stokes equations (run at the 
same friction Reynolds number, and see below for de- 
tails). These modes can be arranged in descending order 
according to their relative energy. Unexpectedly we find 
that the nature of the most relevant modes is unchanged 
in the two cases. On the other hand the energies as- 
sociated with the modes and their relative ordering are 
changed; some modes that are energetic in one flow are 
strongly suppressed (their energy decreases by a factor of 
4) in the other flow, and vice versa. Most importantly, 
the few most energetic modes of the viscoelastic flow con- 
tain a lot more energy that the same number of most 
energetic Newtonian modes. We propose therefore that 
drag reduction should be understood by examining the 
dynamics and relative stability of the energy containing 
modes rather than focusing only on the dissipative end 
of the spectrum, as proposed for example in Q. 

In Sects. H A and B we summarize the FENE-P equa- 
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tions and the numerical approach. In Sect. II C we 



present the essential results regarding the observation of 
drag reduction. In Sect. Ill we review the Karhunen- 



Loeve method for determining the best empirical modes, 
and apply it to the problem at hand. In Sect. ^ we dis- 
cuss the results, demonstrate the invariance of the modes, 
and present the relative ordering. Sect. |v| is devoted to 
a discussion of the findings and of the road ahead. 



II. EQUATIONS OF MOTION AND DIRECT 
NUMERICAL SIMULATIONS 

A. The FENE-P model for dilute polymers 

The addition of a dilute polymer to a Newtonian fluid 
gives rise to an extra stress tensor T(r,t) which affects 
the Navier-Stokes equations: 



du 
~dt 



+ (u ■ V)u 



-Vp + v s V 2 u + V ■ T , 



Vu = 



(1) 



Here u(r,t) is the solenoidal velocity field, p(r,t) is the 
pressure and v B is the viscosity of the neat fluid. The ad- 
ditional stress tensor T is not known exactly, and needs 
to be modeled. There are a number of competing deriva- 
tions, see O, [l^]. In our work we adopt the FENE-P 
model that is derived on the basis of a dumbbell model 
for the polymers jll], [l2| and is known to reproduce the 
phenomenon of drag reduction. In this model ^ is deter- 
mined by the "polymer conformation tensor" R accord- 
ing to 



T(r,t) 



Po 



(2) 



Here 1 is the unit tensor, v p is a viscosity parameter, r p 
is a relaxation time for the polymer conformation ten- 
sor and pq is a parameter which in the derivation of the 
model stands for the rms extension of the polymers in 
equilibrium. The function f(r,t) limits the growth of 
the trace of R to a maximum value p m : 



f(r,t) 



pI 



(3) 



The model is closed by the equation of motion for the 
conformation tensor which reads 



OR, 



■a/3 



Ot 



du a du 7 



- — [f{r, t)R aP - po5 a p] 



(4) 



The model is derived by assuming that the polymer 
can be characterized completely by an end-to-end vector 
distance. Nevertheless the resulting equations could be 
written on the basis of plausible arguments including up 




FIG. 1: Geometry of the channel flow between two par- 
allel planes separated by 2H in the y direction. The mean 
velocity U(y) is oriented along the x axis (streamwise direc- 
tion). The three-dimensional velocity fluctuations are space- 
homogeneous in the x — z plane, where z is usually called the 
span-wise direction. It is customary to use the Fourier de- 
composition in this plane; q = (q x ,q z ) is the corresponding 
two-dimensional wave vector. 



to quadratic terms in gradients and available tensors. For 
our purposes the accuracy of the model in reproducing 
quantitatively all the phenomena of turbulence in vis- 
coelastic fluids is not an issue of concern. We are mainly 
interested in the fact that these equations were simulated 
on the computer in a channel geometry and exhibited the 
phenomenon of drag reduction as discussed below. Our 
aim is to understand drag reduction within the FENE-P 
model. 



B. Direct Numerical Simulations 

A simple flow geometry which exhibits the phe- 
nomenon of drag reduction is channel flow between two 
parallel planes, separated by the distance 2H in the y- 
direction, see Fig. |l|. The computational domain is peri- 
odic in the two directions parallel to the wall (streamwise 
x, spanwise z). The Navier-Stokes equation ([!]) is writ- 
ten in terms of wall- normal component of velocity u y (r, t) 
and the vorticity ui y (r,t), where u> — V x u, 



du y 
~dt 
du 



(U ■ V)Uy 



Of 



~— + V S V 2 Uy + V^y 
(UJ ■ V) Uy + V S V 2 U)y 



(5) 



with 



e Q/ 3 7 the fully antisymmetric Levi-Civita tensor. 

the components of the velocity field 
in the two directions parallel to the x — z plane follow 
from the continuity equation and the definition of ui v , 



Given u y and uj y , 



W u \\ 



dUy 

dy 



Vll X Mil = UJ. 



y> 



(6) 



where the subscript || denotes projection of vectors on 
the x — z plane. The system of Eqs. (||,||) is one of the 
standard formulations used for the direct numerical sim- 
ulation of turbulent channel flows, since the procedure 
yields a solenoidal velocity field to machine accuracy. 
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In light of the periodicity in the x — z plane, it is natu- 
ral to expand the planar components of the velocity field 
in Fourier modes. For the wall-normal direction one uses 
Chebyshev polynomials. The time stepping is carried out 
by a mixed Crank-Nicolson/Runge-Kutta scheme for the 
viscous and the nonlinear terms, respectively. The inte- 
gration in the normal direction is done by Chebyshev tau- 
method and a standard de-aliasing technique is adopted 
for the nonlinear terms. The typical simulations have 
been performed on a computational grid of 96 x 129 x 96 
nodes in a domain of dimensions 2ttH x 2H x 1.2-kH. 
Turbulence is maintained by enforcing the same constant 
pressure gradient for the corresponding Newtonian and 
viscoelastic simulations. 

In discussing the simulations, one notes that the chan- 
nel half-height H is only one of the parameters which 
sets up the external lengthscale. An additional impor- 
tant control parameter is the enforced pressure gradient 
at the wall, which determines the friction. Denoting by 
pointed brackets an average over time the friction param- 
eter is defined as 



, ox 



(7) 



This is the basic control parameter of the flow. By con- 
sidering also the overall kinematic viscosity, Vf = v s + v p , 
the traditional friction Reynolds number is defined as 



Re T = 



u T H 



(8) 



where u T is the friction velocity. As customary in wall 
turbulence, one may introduce the inner, or viscous, 
length scale 



(9) 



The inner velocity scale u T has its counterpart bulk ve- 
locity 



(10) 



Correspondingly we have the outer and inner time scales 
To = H/Uo and Vf/v%. respectively. 

In the sequel all quantities are made dimcnsionless un- 
less stated otherwise. Those made dimensionless with 
respect to the appropriate inner scale will be denoted by 
the superscript +; no special symbol is used to denote 
normalization by an outer scale. 

Finally there are parameters associated with the poly- 
mer. Foremost is the Deborah number which is defined 
as the ratio of the relaxation time r p and a typical time 
scale of the fluid motion, i.e. 



De 



T 



(11) 



The parameter r p = v v /v s measures the relative viscosity 
of the polymers with respect to the Newtonian solvent, 



and the non-linear characteristics of the spring is defined 
in terms of the ratio p^j p\. 

In all our comparisons of Newtonian and viscoelastic 
simulations, the friction Reynolds number is kept con- 
stant, at a typical value Re r = 125. The correspondence 
between the two flows is obtained by fixing the compu- 
tational domain, and choosing for the Newtonian case 
a viscosity equal to the overall viscosity of the solution. 
The parameters chosen for all the viscoelastic simulations 
are De = 25, 7j p = 0.1 and p 2 J p 2 Q = 1000. 



C. Overview of drag reduction 

In this subsection we review briefly the main results 
of the present simulations which demonstrate the phe- 
nomenon of drag reduction. For more detailed descrip- 
tion, see H . The analysis presented in the following sec- 
tions employs the very same DNS. 

In comparing the viscoelastic and the Newtonian flows 
we maintain the friction Reynolds number (^) fixed. We 
reiterate that to achieve this we need to choose the viscos- 
ity of the Newtonian flow properly, since the viscoelastic 
wall stress contains a small viscoelastic contribution, 



dU -=, 
dy 



(12) 



The component of the extra-stress T yx does not vanish 
on the average, contributing in our simulation about 10% 
of the total drag. 

The main observations regarding drag reductions are 
as follows: 

(i) For a fixed mean pressure gradient at the wall, the 
viscoelastic flow exhibits an increased flow rate through 
the channel, see Fig. |^. Shown there are the mean pro- 
files in the streamwise direction, the other means vanish 
by symmetry: 



U x (y) = U(y) = (u x (r,t)), 
U y (y) = U z (y)=0. 



(13) 



The increase in the throughput entails an increased bulk 
Reynolds number, 



Re b = 



UqH 



(14) 



(ii) Some typical length scales increase in the viscoelas- 
tic flow. The data in Fig. ^ can be recast in terms of the 
log-law of the wall, 



U + = -^logy^ 



A, 



(15) 



where k ~ 0.4 is the Von Karman constant. The constant 
A depends on the thickness of the buffer plus viscous sub- 
layers, defined as the distance between the wall and the 
beginning of the log-region. The log-law exists equally 
well for the viscoelastic as for the Newtonian flow with 
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FIG. 2: Mean velocity profiles for the Newtonian and for the 
viscoelastic simulations with Re T = 125. Solid line: Newto- 
nian. Dashed line: Viscoelastic. The straight lines represent 
the classical log-law, Eq. (^). 
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FIG. 3: Two point span- wise correlation of the streamwise 
velocity component, K xx (yo, Z) / K xx (yo,0): viscoelastic flow 
(dashed line), Newtonian flow (solid line). Re T = 125, y$ = 7. 



the same k, but the numerical value of the constant A 
is substantially increased in the former. This thicken- 
ing of the buffer layer is directly related to the increased 
flow rate. The effect of thickening buffer layer can be 
also measured by the increase of the span-wise scale of 
the streamwise velocity fluctuations. Decomposing the 
velocity field into a mean and a fluctuating part 

u a (r,t) = U a (y)+u a (r,t), (16) 

one introduces the correlation tensor 

K a p(r, r') = (u a {r, t)u (r', t)) . (17) 

In Fig. H we show K xx (r,r') for r = (x,y,z), r' = 
(x,y,z + Z). [Note that by homogeneity K xx — 
Kxx{y, Z)]. Fig. | demonstrated the increase in the span- 
wise correlation length of the streamwise velocity fluctu- 
ations which can be defined as 



L X z{y) 



(18) 



(hi) Comparison of the time signals between the vis- 
coelastic and Newtonian cases shows an alteration of the 
characteristic frequencies for the viscoelastic model. This 
corresponds quite well to the experimental data, for ex- 
ample of Luchik et al. (d|, in which a decrease of the 
bursting frequency is observed in drag reducing solutions. 
We do not reproduce these results here. 

(iv) Finally, the root mean square fluctuations change 
significantly as a function of y. Denoting 



U a (y) = ^/(\u a (r,t)\< 



(19) 



we display in Fig. |i| the y dependence of the three compo- 
nents a = x, y, z. The streamwise fluctuations are shown 
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FIG. 4: Velocity fluctuations, normalized by friction velocity 
Ut, for the viscoelastic (dashed line) and the Newtonian flow 



(solid line). In both cases, U : 
dashed and dotted lines, respectively 
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are given by solid 



to increase with respect to the corresponding Newtonian 
flow, while both the span-wise and the wall-normal fluc- 
tuations decrease, in qualitative agreement with available 
experimental results pUfl , 

In conclusion, the FENE-P model is shown to exhibit 
the phenomenon of drag reduction in close similarity 
to experimental observations. In addition, DNS of this 
model provides complete information about the velocity 
field and the covariance tensor field R(r, t) as a function 
of space and time. We thus feel confident that a suf- 
ficiently savvy analysis of this model and its turbulent 
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flow pattern should provide insight into the mechanism 
of drag reduction. 



III. ANALYSIS USING EMPIRICAL MODES 

In this section we present the analysis of the difference 
between viscoelastic and Newtonian flows in term of the 
empirical modes that are sustained in the turbulent flow. 
In choosing a method to extract the modes we are led 
by the desire to find the "best" modes, and "best" in our 
case will mean those that are most energetic. In fluid me- 
chanics the energy is quadratic in the field, so "best" is 
related to closeness in the L 2 norm. As is well known, the 
standard method to find best representation in L 2 norms 
is the Karhuncn-Loeve method. The aim of the Karhuncn 
Loeve method is to provide a set of modes that optimally 
decompose the field of interest, in our case the velocity 
field averaged over time. The approach is guaranteed to 
yield the best set of truncated modes, meaning that the 
field cannot be approximated better (in L 2 norm) by any 
other set of the same number of modes. The method had 
been applied to fluid mechanics in a number of contexts, 
see for example j[4[ [l5|, |l~6| for details and relevant ref- 
erences. We first adapt the Karhuncn-Loeve method to 
the present context, and then discuss the results of the 
analysis. 

The core object for the analysis is the simultaneous cor- 
relation function of the velocity fluctuations which was 
introduced already in Eq.(|l7j). Due to stationarity this 
object is time independent, and due to homogeneity in 
the x and z direction, we can write 



K a = K a p(x' - x, z' - z; y, y') 



(20) 



In the translationally invariant directions x—z one cannot 
do better than Fourier decomposition. Accordingly we 
consider the partly decomposed object Q(q\y, y') defined 
as 

Q(l\y,y') = -^-^K aj3 (x - x : z' - z\y,y') (21) 

II P|| 

X exp{-i[q x (x' -x) + q z (z' - z)]} . 

We denoted the discrete two dimensional wave-vector in 
the x - z plane as q = {q x ,q z ) = (q x 2n/A x , q z 2n/A z ) 
where q x ,q z are integers. _/V|| = N X N Z and the sum is 
taken over the discrete set of x — z points with rn = 
(jxA x /N x ,j z A z /N z ) the x - z projection of r. 

The non-trivial empirical modes are obtained from 
the remaining dependence on for each given pla- 

nar q. The Karhuncn-Loeve method consists of finding 
the eigenfunctions ^ a p(q,p\y) which solve the eigenvalue 
equation 



Q a p{q\y,y')^fi(qMy'W = E(q,p)^f a (q,p\y) 



Here E(q,p) is the energy associated with the mode 
(q,p), ordered in decreasing energy order, with a given 
q wave- vector in the wall parallel directions, and p mode 
number in the y-direction, orthogonal to the plane walls. 
Denoting by £ the total energy in the system, we have 



<J,P 



(23) 



The modes labeled by (q, p) can be relabeled in order of 
decreasing energy by introducing an index n. In this no- 
tation the mode whose associated energy is E n = E(q, p) 
are denoted by 



V a (n\y) = V a (q n ,Pn\y), 



(24) 



where (q n ,Pn) is the label of the n th mode. It is useful 
to introduce also a set of fields indexed by n, 



$ a {n\r) = ^ a {n\y) exp [i(q n ■ r\\)] 
These function arc orthogonal in the sense 



1 f 1 

— 22 / <fo$ a (n|i/,r||)$£(n'|y,ry) 

II -r-,, •* — 1 



= <5n 



(25) 



(26) 



The fluctuation velocity field can be now expanded in 
terms of these fields, 



t (t,r) = ^2a n (t)$ a (n\r). 



(27) 



The main advantage of the Karhunen-Loeve method is 
that any finite truncation of this expansion can be shown 
to yield a best approximation for the velocity field (in the 
L 2 norm). Moreover, due to the symmetry of the corre- 
lation matrix the modes are orthonormal (after suitable 
normalization), and correlation functions in the basis of 
these modes are diagonal: 



(a„a*,) — E n 5 n 



(28) 



Similarly one can have an optimal representation of the 
correlation matrix itself as 

K a p = ^2E n ^ a (n\y)^p(n\y')exp iq n ■ (r\\ -rj) 

n 

(29) 

In analyzing the DNS the modes arc determined by 
using an expansion in terms of Chebyshev polynomials 
Tfc(y) in the y direction. Denote yj the j" 1 node in the 
list of N y Chebyshev nodes, 



and express 



cos j 



j = 1, • 



(22) 



*a(q,p\y) = * a (q,p\k)T k (y) 



(30) 



(31) 



k=0 
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Then Eq. 



is discretized as 



IV. ANALYSIS AND RESULTS 



Ny-l 

E 

k=0 



A af} (q\j,k)y/3(q,p\k) 



Ny-l 

E(q,p) **(q,p\k')T k ,( yj ) . (32) 

k'=0 



where 



A a [}(q\j,k) 



N v 

E Qccp{Q\VhV3') T k{Vj')Wj' 



(33) 



with Wj the integration weights. 

For a given q, Eq. (^) results in a linear algebraic 
eigenvalue problem of order (N y x N y ) , with the p th eigen- 
vector given by ^ a (q, p\k), k — 0, . . . , N y and the corre- 
sponding eigenvalue E{q 1 p). Since we need to solve for 
every discrete wave-vector q, there exist in total N X N Z 
eigenvalue problems to be solved in order to determine 
the full set of eigenfunctions. From the discrete eigenvec- 
tors, the corresponding discrete modes are constructed 
according to the discrete analog of Eq. (E5h , 



® a (q,p\r\\,yj) = * Q (g,p|yj)exp [i(q n • ry)] , (34) 

where Eq. (|3l| ) is used to evaluate ^f a at the Chebyshev 
nodes. The discrete modes are normalized according to 
the discrete version of Eq. ( p6| ) , 



II j=l r y 

In practice we ran DNS with and without polymers 
for 5000 large-eddy turnover times T , see Sec. II B, in 
statistically stationary conditions. The time step used 
for t ime- advancement, in terms of viscous units Vf/v% 
Sec. ~ 
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is Dt + = 0.05. Since the same pressure gradi- 
ent is enforced in the two cases, the friction velocity is 
identical, while the bulk velocity increases by 24% in the 
viscoelastic simulation. This corresponds to a reduction 
of the large-eddy turnover time, To. 

We have collected iVfieids fields, iVfieids — 100, displaced 
in time by 50 T . These fields were used to construct the 
discrete correlation function, Eq. (Jl7|), at the nodes of 
the computational grid. Eigenvalues and eigenvectors of 
its discrete Fourier transform in the wall parallel direc- 
tions has then been evaluated, according to Eq. (^) . The 
modes corresponding to the computed eigenvectors have 
been arranged in order of decreasing eigenvalue. From 
Eq. (^||), the eigenvalues correspond to the average en- 
ergies in the considered modes, so that the chosen ar- 
rangement is in fact an ordering in terms of decreasing 
energy content. When needed, the index of the energetic 
ordering of the Newtonian and viscoelastic simulations 
will be denoted by n N and n VE respectively, dropping 
the subscript when n refers to both cases. 



A. The dominant modes are approximately 
invariant 

The first discovery in the study of the empirical modes 
was admittedly surprising for the present authors, and in 
hindsight very serendipitous for the discussion of drag re- 
duction. We found that the dominant modes are approx- 
imately invariant. In other words, the same modes that 
carry a sizable fraction of the energy of the viscoelastic 
flow appear essentially unchanged in the Newtonian flow, 
having practically the same spatial y-dependence. We 
will first demonstrate this surprising finding, and later 
focus on the difference between the two flows. 

To discuss meaningfully the correspondence between 
the empirical modes in the two cases we seek a criterion 
of matching the viscoelastic modes to the correspond- 
ing Newtonian modes. This can be done easily, since 
both sets of modes form a complete orthonormal basis 
for solenoidal fields in the same domain. Each viscoelas- 
tic mode can then be expanded in terms of the Newtonian 
set, 

C E (g ) Pv E |y)=E A (9'Pv E IPN)*a(9.P N |2/) • (36) 

Pn 

The complex amplitude A( q, p VE |p N ) is given by the pro- 
jection of the viscoelastic mode p VE on Newtonian mode 
p N , with identical wall-parallel wave- vector, q: 



^(<?,PveIPn) 



'VE I 



(37) 



2N, 



]T*H<7, P VE |2/,)*f(g, pJy.H 



We find that all the most energetic modes of the vis- 
coelastic flow have one amplitude whose magnitude is 
close to unity. For example we plot the absolute magni- 
tude |A(q,p VE |p N )| as a function of p N in Fig. || for the 
first three most energetic viscoelastic modes. Each of 
these modes displays an amplitude maximum very close 
to unity, implying that it matches well a single Newto- 
nian mode. This procedure furnishes a correspondence 
between viscoelastic and Newtonian modes: a given vis- 
coelastic mode corresponds to the Newtonian mode for 
which the amplitude |A(<7,p VE |p N )| is maximal. This un- 
ambiguously associates a single Newtonian mode to each 
viscoelastic mode. The value of the maximal amplitude, 
hereafter called matching parameter, gives a quantitative 
estimate of the difference in shape between two corre- 
sponding modes: A matching parameter equal to one im- 
plies absolute identity between the two modes in the en- 
ergy norm. The correspondence in the spatial structure 
can then be verified by direct inspection. For instance, 
matching modes for the cases discussed in Fig. @, where 
the matching parameter is well above 0.9, are almost in- 
distinguishable. Physically, the matching parameter rep- 
resents the fraction of the energy in the viscoelastic mode 
that is ascribed to the matching Newtonian mode. 
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FIG. 5: Projections of the three most energetic viscoelas- 
tic modes on the basis of the Newtonian best modes with 
the same q, Re T = 125. Panel a: the most energetic Vis- 
coelastic mode, q = (0,3), p = 1. This modes fits very well 
the most energetic mode of the Newtonian flow which is the 
same, (0,3,1). Panel b: the second most energetic Viscoelastic 
mode, q — (0, 2), p = 1. This modes fits very well the sixth 
most energetic mode of the Newtonian flow (0,2,2). Panel c: 
the third most energetic Viscoelastic mode, q = (0, 1), p = 1. 
This mode fits well the fourth most energetic Newtonian mode 
(0,1,1) 



Fig. |^ shows the matching parameter for the first thirty 
most energetic viscoelastic modes (n VE < 30). All the 
VE-modes except the 23rd and 29th have a matching 
parameter above 0.9. We conclude that at least as far 
as the most energetic modes are concerned, the spatial 
structure of the modes is almost not altered by the poly- 
mer. Loosely speaking this can be expressed by saying 
that the modes are approximately invariant, i.e. fixed 
in shape. For practical purposes the difference between 
the Newtonian and viscoelastic empirical modes can be 
safely disregarded. Next we discuss what changes from 
Newtonian to viscoelastic turbulence. 
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FIG. 6: Matching parameter of the first 30 most energetic vis- 
coelastic modes with corresponding Newtonian modes, Re r = 
125. 




FIG. 7: Iso-surface of the streamwise component of the first 
most energetic mode, $i(n = l|r) which is the same for the 
Newtonian and the viscoelastic flows, i.e. the mode q = (0, 3), 
p = 1. In both cases, Re T = 125. Note that this mode is a 
non-propagating roll mode. 



B. The energy and the relative ordering change 

In light of the first discovery the second may be al- 
ready anticipated: although the dominant modes hardly 
change, their energies and relative energy ordering 
change very significantly. We propose that understanding 
the change of energies and relative ordering of the modes 
will take us a long way in understanding drag reduction. 
We begin by comparing the most dominant modes in the 
two respective flows. Each empirical mode is identified by 
three numbers, (q x , q z ), and p, where as explained above, 
p corresponds to the index of the energetic ordering for 
fixed wave- vector q. 



1. The most dominant mode 

In Fig. we present pictorially the mode (0,3,1) which 
is the most energetic for both the Newtonian and vis- 
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FIG. 8: Portrait of the most energetic Newtonian mode, 
(0,3,1): Upper panel - real part of the fluctuating velocity 
profile *P(l|y), lower panel - imaginary part of vE'(ljy). Solid 
lines - ^ x (l\y), dashed lines - ty y (l\y) and dash-dotted lines 

coelastic simulation. Its amplitude a±(t), 
a i(t) = u a (t,yj, r||)$*(n= %•, r\\)wj , 

v j=l r y 

is real in both flows, implying that the mode does not 
propagate, neither in the x nor in the z direction. More- 
over, for this particular mode q x = 0, i.e. the field 
$ Q (n = l|r) is constant in x. It belongs to the class 
of non- propagating roll- modes discussed in pl|. The 
figure presents an iso-surface of the streamwise com- 
ponent of the velocity field associated with the mode, 
$> x (n = l|r) = const. The span- wise wavelength of the 
mode is X z = A z /q z = 2itH/3, and it appears to be con- 
fined relatively close to the channel walls. 

A further pictorial presentation of the Newtonian mode 
(0,3,1) is shown in Fig. ||, where the real part of the mode 
(top panel) and imaginary part (bottom panel) are plot- 
ted separately. From the plots it is apparent that the 
mode is more or less localized in the vicinity of the walls, 
as already commented. The phase difference between 
the various components is also worth mentioning. Com- 
paring the top and the middle panels, the stream-wise 
and span-wise components, ^ x and <F z are real for all 
?/, while the wall-normal component, 'Fy, has a constant 
phase lag of ir/2 with respect to the other two, i.e. it 
is purely imaginary. We note that although the first, 
most dominant mode, is the same for the two flows, the 
actual energy associated with this mode is twice larger 
in the viscoelastic mode. This is typical for all the lead- 
ing modes in the viscoelastic flow as compared with the 
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FIG. 9: Portrait of a typical viscoelastic high order mode, 
(0, 3, 15). Notations as in Fig. |. 

Newtonian flow; this central point of distinction between 
the two flows will be addressed in the next sub-section. 
The agreement between the leading modes ends with the 
first one; the majority of higher modes change in relative 
position in the energy descending ordering, as explained 
below. Note that all the seven leading modes in both 
flows are associated with q x = 0. This is surely due to 
the relatively short channel length in our computation. It 
is known that in longer channels the "roll" modes become 
oscillatory modes with a finite value of q x . 

For the sake of completeness we present in Fig. ^ 
the portrait of a higher order viscoelastic mode, namely 
q = (0, 3), p = 15. We see that for this mode 'F z is 
essentially purely imaginary, whereas ^ x and ^ y are es- 
sentially purely real. In fact, this particular mode is very 
low in energy, and we present it just to demonstrate that 
the numerics is still not noisy even for rather low energy 
modes. 



2. Sub-dominant modes and energy redistribution 

The full comparison between the first 30 modes of the 
Newtonian and viscoelastic flows can be seen in Table 1, 
in which these modes are listed together with their energy 
(in percent of the total sum of energies). Also included 
in the table is the cumulative sum of energies up to the 
mode listed. 

To understand the Table we recall that for each value of 
q the modes are ordered by their energy and labeled by p. 
The over-all energy label is n N or n VB . In the columns of 
the viscoelastic modes, instead of writing n VE in obvious 
increasing order, we provided the value of the Newtonian 
index n N for which the viscoelastic mode has the high- 
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TABLE I: Energy content of the first 30 eigenfunctions of the 
channel flow with and without polymers. Re T = 125. 

est matching parameter. Thus for example the second 
leading viscoelastic mode with n VE = 2 matches almost 
perfectly the sixth Newtonian mode n N = 6, etc. One 
glaring difference between the two flows that stands out 
from Table 1 is the energy concentration in the few most 
dominant modes of the viscoelastic flow compared to the 
flat distribution of energy in the Newtonian flow. For ex- 
ample, the first, second and third dominant viscoelastic 
modes have all twice the energy of the corresponding first, 
second and third dominant Newtonian modes. The sum 
of the first four viscoelastic energies contain as much as 
as the first ten Newtonian modes. The first 15 viscoelas- 
tic modes already contain more than 50% of the energy, 
whereas one needs to collect as many as 80 Newtonian 
modes to reach the same fraction of the total energy. In 
Fig. [l(] we display the sum £{n) = X^=i as a func- 
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FIG. 10: Energy sum £(n) = XT=i f° r ^ ne Newtonian 
(solid lines) and for the visco-elastic (dashed lines) flow. The 
sums are computed in their respective energetic ordering. The 
Reynolds number is Re T = 125. 



tion of n for both flows, where the sum is computed in the 
respective energy ordering. We note that the difference 
between the two curves is established within the first ten 
modes; for larger values of n the lines are almost parallel, 
indicating a similar energy distribution between the less 
dominant modes. 



C. The energy spectrum 

An interesting and illuminating way to discuss the dif- 
ference between the two flows is provided by the energy 
spectrum. It was proposed in [|] that the main difference 
between the spectra of the two flows should appear in the 
position of the dissipative scale that separates a spectral 
power-law from exponential decay. The increase of the 
dissipative scale was indeed observed in recent DNS of 
homogoneous isotropic turbulence in the FENE-P model 
JT7| . However, the present results indicate that impor- 
tant changes involve the energy containing scales. 

In Fig. |ll| the energy is plotted for the two flows as 
a function of n. The Newtonian case displays a spectral 
plateau for the most dominant modes which crosses over 
to a power law for n N > 8. In the case of the viscoelastic 
flow the plateau is dramatically higher, but also the power 
law changes its slope. It appears that the whole spectral 
curve is tilted in favor of the energy containing modes 
and on the expense of the lower energy modes and the 
dissipative modes. 

The difference between the power laws can be made 
clearer by comparing with what can be obtained from the 
Kolmogorov law for high fc-vectors. For relatively large k 
we can expect the flow to isotropize |l|, |2(| , and Fourier 
modes would again become "best". In this asymptotic 
situation (neglecting intermittency corrections) the spec- 
trum E{k) is expected to be E{k) oc fc -11 ' 3 . In an 
isotropic environment we also expect that n is propor- 
tional to the volume of the sphere of radius k, i.e. n oc k 3 . 
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FIG. 11: Log-log plots of the energy distribution E n vs mode 
index n for the Newtonian case (solid line) and the viscoelastic 
case (dashed line). The rough K41 prediction (^) is shown 
by the dashed-dotted line. The Newtonian and visco-elastic 
dependencies are shown in their own energetic ordering. The 
Reynolds number is Re T = 125. 



Thus we can expect for large values of n 
E n oc n -11 / 9 , n large . 



(39) 



The law proposed in Eq. (|3^) is displayed as the dashed- 
dotted line in Fig. [ll[ We see that it is in rough agree- 
ment with the power law section of the Newtonian spec- 
trum, but it is certainly not in agreement with the vis- 
coelastic spectrum which, as said before, is becoming 
steeper on the whole. This is a clear demonstration of 
the increase in energy of the energy containing modes 
on the expense of the others. We propose that even the 
power law section of the viscoelastic spectrum will not 
have a "universal" slope, but rather a slope that depends 
on the concentration of the polymer and the degree of 
drag reduction. 

What emerges from this study is that understanding 
drag reduction lies in the reordering and energy redistri- 
bution of the energy containing modes. To examine these 
phenomena further we consider now the energy contents 
of the viscoelastic modes ordered by n N instead of their 
own ordering. This plot is a very vivid graphic repre- 
sentation of the data of Table 1, stressing very strongly 
the concentration of energy in the dominant viscoelas- 
tic modes, which are however rearranged in dominance 
compared to the Newtonian case. Next we want to reit- 
erate that one should focus on the most energetic modes. 
Noticing from Table 1 that all the dominant modes in 
both flows are associated with q x = 0, we consider next 
these modes as a function of q z for p = 1, 2. In Fig. [l^ we 
show the energy of these modes. We note that the dra- 
matic redistribution of energy occurs only for the most 
energetic modes with p — 1; The less energetic modes 
are not affected much. Already the modes with p ~ 2 are 
seen in the figure to be affected in a negligible way. In our 
opinion this is a clear message that to understand drag 
reduction we need to understand the rearrangement of 
the energy containing modes. For our flow configuration 
the most relevant are the modes which are space homo- 




FIG. 12: Plot of the energy content (in percentage of the 
total energy) for the Newtonian (solid line) and viscoelastic 
(dashed line) empirical modes. Both cases are plotted as a 
function of n„, . 




FIG. 13: Energy of the most energetic modes for a given 
q x = as a function of q z . p = 1 modes are represented by 
circles and p = 2 modes by triangles. They are connected 
by solid lines for Newtonian modes, and by dashed line for 
viscoelastic modes. 



geneous in the spanwise direction. It should be stressed 
however that different geometries, and even channel ge- 
ometry with different aspect ratios, may bring forth other 
modes as the most relevant ones. Nevertheless we ex- 
pect that drag reduction would always be associated with 
a substantial increase in the energy containing modes, 
whatever these are for a given flow configuration. 



V. CONCLUSIONS 

In this paper we initiated a systematic study of drag 
reduction on the basis of DNS of the FENE-P model. 
We investigated simulations of Newtonian and viscoelas- 
tic flows in channel configuration at the same friction 
Reynolds number. Our main aim is to understand the 
mechanism of drag reduction. Since drag reduction in- 
volved modifications of the mean flow and of the large 
scale gradients in the flow, we are motivated to under- 
stand more the energy containing modes, rather than fo- 
cusing only on phenomena of small scales. To this aim we 
have found first the list of empirical modes that represent 
the velocity field in an optimal way in an energy decreas- 
ing ordering. The first important discovery was that this 
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list contains the very same modes for the Newtonian and 
viscoelastic flows. We propose that this finding will offer 
a huge simplification in any future theory of drag reduc- 
tion. While we cannot offer a clean explanation of this 
finding, we can proceed at this time taking the approxi- 
mate invariance of the modes as an empirical fact. What 
needs to be understood is just the energy distribution 
and reordering of the modes in the viscoelastic case. 

We should stress that the point of view proposed here 
differs in a fundamental way from the approach presented 
for example in Refs. Q |l8). The thinking there fo- 
cuses on the energy cascade in the turbulent flow, and on 
the modification of the small dissipative scales. Roughly 
speaking, the maximal gradients of the velocity are es- 
timated from balancing the RHS of Eq. (^). Neglect- 
ing the statistical correlation between the conformation 
tensor and the velocity gradient, one can estimate the 
maximal velocity gradient by (du/dr) ~ 1/t p . Thus the 
matching of the polymer relaxation time with an eddy 
turn over time is used to predict a decrease in the max- 
imal velocity gradient which is interpreted as drag re- 
duction. We take an exception to this approach. First, 
a careful analysis of the space dependence of the con- 
formation tensor shows that it is highly correlated with 
the velocity gradients (see for example ||, [ij], ^l| . Thus 
the estimate taken above is questionable at best. But 
moreover, we have shown that the main changes between 
the Newtonian and the viscoelastic modes occur in the 
energy containing modes. The energy containing modes 
are highly anisotropic, they are not Fourier modes, and 



their connection to the small scales where the isotropised 
Kolmogorov picture is tenable is very unclear. A theory 
that assumes a K41 spectrum down to a modified viscous 
scale does not appear tenable for the FENE-P flow, as 
seen in Fig. O. We propose that a theory of drag re- 
duction entails an understanding of the relative energy 
of the modes that characterize the largest scales of the 
flow. 

To understand the relative ordering and the relative 
energy of the most dominant modes one needs to study 
the energy intake by these modes from the mean flow 
H 0, the energy exchange between the modes, and the 
energy exchange with the viscoelastic subsystem repre- 
sented by the conformation tensor field R(r,t) Jl7], pi| . 
Such an investigation calls for measuring additional sta- 
tistical objects like 3rd order correlation functions. The 
necessary simulation are in progress and the results will 
be published elsewhere. 
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